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Abstract 

We report a Monte Carlo simulation of the 2D Edwards- Anderson spin glass model 
within the recently introduced multicanonical ensemble. Replica on lattices of size 
up to L = 48 are investigated. Once a true groundstate is found, we are able to give a 
lower bound on the number of statistically independent groundstates sampled. Tem- 
perature dependence of the energy, entropy and other quantities of interest are easily 
calculable. In particular we report the groundstate results. Computations involving 
the spin glass order parameter are more tedious. Our data indicate that the large L in- 
crease of the ergodicity time is reduced to an approximately power law. Altogether 
the results suggest that the multicanonical ensemble improves the situation of simula- 
tions for spin glasses and other systems which have to cope with similar problems of 
conflicting constraints. 
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1 Introduction 



The theoretical understanding of spin glasses (for reviews see |]T], has remained a great 
challenge. In particular the low temperature limit leaves many open questions about the ef- 
fects of disorder and frustration. For instance, it has remained controversial whether Parisi's 
[0] mean field theory provides the appropriate description for 3D spin glasses. The attrac- 
tive alternative is the droplet model [Q, which in turn is equivalent to a one parameter 
scaling picture . The simplest spin glass system to study such questions numerically is the 
Edwards-Anderson model. In its Ising version it is described by the Hamiltonian 

H = -Y. JijSi^v (1) 

<ij> 

where the sum goes over nearest neighbors and the exchange interactions Jij = ±1 between 
the spins Si = ±1 are quenched random variables. In our investigation we impose the 
constraint J2 Jij = for each realization. Despite its simplicity the model is supposed to 
be sufficiently realistic to catch the physics essence correctly. Recent simulations of the 
3D model in a magnetic field support the mean field picture. However, one may well argue 
that sufficiently low temperatures on sufficiently large systems have not been reached. For 
previous simulations of the Edwards- Anderson model without magnetic field in 2D, 3D and 
4D, see 

Low temperature simulations of spin glasses suffer from a slowing down which is likely 
to increase with a high power law or even exponentially fast with lattice size. The reason 
is that one has to sample many independent states separated by energy barriers which may 
grow with lattice size. To illustrate the problem, let us consider a simple ferromagnet: the 
2D Ising model on a 50 x 50 lattice. In Figure 1 we give its magnetic probability density. 
The two distinct branches below the Curie temperature are associated with free energy 
valleys in configuration space, each of which defines a (pure) thermodynamic state. The 
notations phases or ergodic components are also used. At temperatures below the Curie 
point the ergodicity time^ r| increases exponentially fast with lattice size, asymptotically 
like exp[/'^(/3)L^~^], where is the surface free energy. Therefore, on large lattices at 
sufficiently low temperature the simulation of the system will, given a reasonable finite 
amount of computer time, never tunnel from one phase to the other. Besides for particular 
problems, like the determination of the order-order surface free energy |T^, |n[], this lack 



of tunneling does not impose a major handicap on Ising model simulations. The reason is 
that the two configuration space valleys are related by the exact symmetry Sj —Si of the 
Hamiltonian. Exploring one valley by means of a simulation yields also all the properties of 
the other one and, hence, allows to overlook the entire system. 

The situation is much more involved for spin glasses. For low enough temperature the 
system is supposed to split off into many thermodynamic states, separated by similar tun- 
neling barriers as the two pure states of the Ising model. However, unlike in the Ising model 
the states are not related to each other by a symmetry of the Hamiltonian. Rather they 
appear because of accidental degeneracy which in turn occurs because of randomness and 
frustration of the system. For computer simulations this means that one would like to explore 



"'^As will become clear in the next section, it is for the present investigation of spin glasses more appropriate 
to use the term ergodicity time r*^, instead of tunneling time r* which is appropriate in the context of surface 
tension investigations. 
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Figure 1: Ising model magnetic probability density from a 50 x 50 lattice. 



many independent configuration space valleys while keeping track of their relative weights. 
The groundstate energies associated with these valleys may or may not be degenerate, but 
it should be noted that even if they are not degenerate, tunneling between the valleys would 
still be governed by the energy barriers. The physics of these barriers is far less well un- 
derstood as in the ferromagnetic case. As detailed finite size scaling (FSS) studies do not 
exist, it is unclear to us to what extent these barriers depend on the system size, whereas 
the temperature dependence has been investigated We use the notation bifurcation tem- 
perature (bifurcation point) for the temperature at which the spin glass configuration space 
(phase transition or not) begins to split off into a number of valleys which are well separated 
by energy barriers. In the present paper we suggest that the increase of r|; can be reduced 
to a fairly decent power law by performing a simulation which covers in a single ensemble a 
whole temperature range from well above to far below the bifurcation point. The appropriate 
formulation is provided by a generalization of the multicanonical ensemble which we 
introduce in section 2. To test the approach we have performed multicanonical simulations 
of the 2D Ising ferromagnet and then of the 2D Edward- Anderson Ising spin glass model, and 
our numerical results are reported in section 3. We concentrate on ground state properties 
what is a kind of worst case scenario for the performance of the multicanonical simulation. It 
should be noted that Figure 1 does not exploit exact results like the symmetry Si —Si, but 
relies on a multicanonical simulation, what explains the slight asymmetry between the two 
branches. In section 4 the multicanonical performance in comparison with other simulation 
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methods is evaluated in some detail. Our conclusions are summarized in section 5. A short 
account of the present work has already been reported [|l^ . 



2 The Multicanonical Ensemble 

Ever since the pioneering paper by Metropolis et al.[14| most MC simulations concentrated 
on importance sampling for the canonical Gibbs ensemble. It has always been well-known, 
for instance [jl5[, that it is allowed to choose phase-space points according to any other prob- 



ability distribution, if it is convenient. However, a systematic reasoning for a better than 
canonical choice has rarely been put forward.^ It is our suggestion that in a large class of 
situations, in particular those where canonical simulations face severe ergodicity problems, it 
is more efficient to reconstruct the Gibbs ensemble from a simulation of a multicanonical en- 



semble than simulating it directly. In canonical simulations configurations are weighted 
with the Boltzmann factor 

PBiE) = expi-pE). (2) 

Here E is the energy of the system under consideration, and in this paper we use the notation 
/3 in connection with the canonical ensemble. The resulting canonical probability density is 

P,{E) ~ n{E)PB{E), (3) 

where n{E) is the spectral density. In order of increasing severity problems with canonical 
spin glass simulations are: 

• i) Simulations at many temperatures are needed to get an overview of the system. 

• ii) The normalization in equation (3) is lost. It is tedious to calculate important 
physical quantities like the free energy and the entropy. 

• in) The low temperature ergodicity time r| diverges fast with lattice size (either ex- 
ponentially or with a high power law). The relative weights of pure states can only be 
estimated for small systems. 

Let us choose an energy range £^min < E < E^ax and define for a given function P{E) 
the function a{E) by the recursion relation (with the Hamiltonian (1) the energy changes in 
steps of 4) 

a{E - 4) = a{E) + [PiE - 4) - P{E)] E, a{E^^^) = 0. (4) 

The purpose of the function a{E) is to give (3{E)~^ the interpretation of an effective tem- 
perature. The multicanonical ensemble |[T^ is then defined by weight factors 

Pm{E) = exp [-/5(E)E + a(E)] , (5) 

where P{E) is determined such that for the chosen energy range the resulting multicanonical 
probability density is approximately fiat: 

Pmu{E) = Cmu n{E)PM{E) ^ const. (6) 



^Notable may be microcanonical simulations, but for the ergodicity problems on which we focus here they 
perform even worse than the canonical approach does. 
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In the present study we take -Emax = {P{E) = for E > -Emax) and -Emin = the ground 
state energy of the considered spin glass reahzation. 

In contrast to first order phase transitions |]T2|] , where finite size scaling allows an accurate 



determination of the multicanonical parameters from previously simulated smaller systems, 
the situation is now more involved. Already changing the realization, i.e., the quenched 
random variables Jij, for fixed lattice size leads to a situation which requires to calculate 
new multicanonical parameters from scratch. Our main experience with such calculations 
is that they can be done and that they are less problematic than one might superficially 
expect. For instance, a multicanonical function l3{E) can be obtained via recursive MC 
calculations. One performs simulations with (3'^{E), n = 0,1,2,..., which yield probability 
densities P^{E) with medians -^median- ^ < -^min < ^median the probability density 
P^{E) becomes unreliable due to insufficient statistics, caused by the exponentially fast fall- 
off for decreasing E. We start off with n = and (3^{E) = 0. The recursion from n to n + 1 
reads 

^E) for E > 



+0.25X ln[P"(E + 4)/P"(E)] 

for E median > E > E^,-^; 



Here the simulation may be constrained to E < by rejecting all proposals with 

energy E > -Emldian' but one has to be careful with such bounds in order to maintain 
ergodicity. The recursion is stopped for m with = being groundstate. 

Starting with this simple approach we have explored several more sophisticated variants. 
Considerable speed-ups and gains in stability could be achieved. The CPU time spent to es- 
timate the multicanonical parameters was 10% to 30% of the CPU time spent for simulations 
with the final set. Nevertheless, our determinations of the multicanonical parameters have 
remained kind of unsystematic. We are not yet able to report a theoretically sound optimized 
automatic procedure, although we are convinced that this will be the final outcome. 

Once the functions I3{E) and a{E) are fixed, the multicanonical simulation exhibits a 
number of desirable features: 



i) By reweighting |T^, with exp[—f3E + (3{E)E — a{E)] the canonical expectation 
values 



where 



O0) = Z{P)-'J20{E)n{E) exp(-/3E), (8) 



Z0) = Y^niE) exp{-PE) (9) 

E 



is the partition function, can be reconstructed for all /3 in a range Pram < /3 < /3max, 
where Pram = PiEraax) and Praax = /^(-E'min) follow from the requirement -Emax > E{P) > 
Emin, and E{P) follows from (8) with 0{E) = E. This feature inspired the name 
multicanonical ensemble. With our choice -Emax = and -Emin = E^ groundstate, 

Pmin = and /3max = CO folloWS. 
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ii) The normalization constant Cmu in equation (6) follows from Z{0) = J2e n{E) = 2^ , 
where is the total number of spin variables. This gives the spectral density and allows 
to calculate the free energy as well as the entropy. (Remember, /3 = is included in 
our choice of the multicanonical ensemble.) 

Hi) We conjecture that the slowing down of canonical low temperature spin glass sim- 
ulations becomes greatly reduced. For the multicanonical ensemble it can be argued 



12| that single spin updates cause a ID random walk behavior of the energy E. As 



-f'max — -f'min ~ V, ouc uccds V"^ Updating steps to cover the entire ensemble. For first 



order phase transition the observed slowing down [|l^, |TT[ was only slightly worse than 
this optimal behavior. Our present MC data show more drastic modifications for spin 
glass simulations. See section 3 for further details. 



To quantify our discussion of the slowing down, we have to define the ergodicity time 
r£. Roughly speaking it is the CPU time needed to collect independent configurations, 
here for groundstates which are our main interest. Regarding the definition of r|, an L- 
independent over-all factor is free. For the purposes of this paper we define r£ as the 
average number of sweeps needed to move the energy from -Emax to -Emin and back. A 
sweep is defined by updating each spin on the lattice once (in the average once if random 
updating is implemented). We denote by Ur the number of "tunneling" events with respect 
to the ergodicity time (suppressing now the obvious L subscript). Assuming that the correct 
groundstate is found, Ur gives a lower bound on the number of independent groundstates 
sampled. This follows from another trivial, but remarkable property of the multicanonical 
ensemble: Each time a sweep is spent at I3{E) = 0, the memory of the previous Markov chain 
is lost entirely, and a truly independent new series of configurations follows. The condition 
E > -Emax is appropriate to substitute for the somewhat too strict constraint of an entire 
sweep at /5 = 0. As a corollary: with a disordered starting configuration the multicanonical 
ensemble is immediately in equilibrium. 

The energy density e(/3), the specific heat c(/3), the free energy density /(/3) and the 
entropy per spin s{/3) follow in a straightforward manner by reconstructing the canonical 
ensemble (8,9). In this paper we go for the extreme and concentrate on the zero temperature 
(/3 oo) limit. Let N = L x L he the total number of spins (A^ = in our notation). The 
groundstate energy density is e° = E'^/N, and we obtain its entropy per spin from (8) as 

= S'^/N = ln[n{E")]/N. More complicated are calculations which aim at the spin glass 
order parameter q. One way to define g is as the overlap of two statistically independent 
replica 

1 ^ 

= M^'l^- (10) 

i=l 

Here we consider equilibrium configurations with respect to the canonical ensemble, sj are 
the spins corresponding to configurations of replica one and sf are those corresponding to 
replica two. Different replica have identical quenched random variables Jij and in practice 
the independence is achieved by creating uncorrelated disordered starting configurations. 
The problem is now that we slow down by at least another factor 1/V when we simply 
simulate each replica with respect to the multicanonical ensemble. The reason is that the 
configurations sj and will normally be at vastly different effective temperatures l3{E^)~^ 
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and P{E^)~^, with then being suppressed by the reweighting procedure at all tempera- 
tures /3. Here and E"^ denote the energy replica one and two, respectively. To by-pass the 
problem, we decided to constrain the multicanonical simulation of the second replica to an 
appropriate energy range E^ — AE < E"^ < E^ + AE and to update until E"^ hits the value 
E^. In this way, we measure a microscopic spin glass order parameter q{E), which can be 
reweighted without problems. In the infinite volume limit q{E{(3)). The spin glass 

susceptibility density Xq ^^nd the Binder parameter are then defined as usual 



Xq = (g^) and Bg = ^ 



3- ^^'^ 



(11) 



Again, we are mainly interested in the zero temperature results, denoted by Xq and B^. The 
constrained MC of the second replica eats the bulk of the CPU time. In the average one 
needs several updating sweeps to find E"^ = E^, because one has to choose AE generous 
in order to avoid ergodicity problems. Although our procedure works, and our results of 
section 3 are based on it, we do not recommend it for future use. The better approach 
seems to be to use instead of (10) an overlap function which is defined with respect to the 
groundstate ensemble [||]. As tunneling with respect to our ergodicity time separates truly 
independent groundstates, one can even avoid the second replica entirely, while keeping the 
physics goals in essence unchanged. Unfortunately, we have not yet practical experience with 
this approach. 

To conclude this section let us comment on the connection with simulated annealing 
[0. A major shortcoming of simulated annealing is that the connection with canonical 
equilibrium configurations gets lost. The multicanonical ensemble overcomes this problem. 
The price paid is random walk-like energy changes in contrast to more directed changes 
in case of simulated annealing. Consequently, the multicanonical ensemble is confined to 
more limited temperature gradients than simulated annealing. Whether this is really a 
disadvantage remains to be explored, as the smaller temperature gradients will increase the 
ability of the simulation to find its way around (or in and out of) metastable states. In 
simulated annealing the temperature gradients are free parameters, whereas in our case they 
are fixed by the ensemble. 



3 Simulations 

All our numerical calculations were performed on the SCRI cluster of RISC workstations. 
Aiming mainly at simplicity, we developed a simple Metropolis type program. Different 
realizations of the spin glass system correspond to statistically independent sets of quenched 
random variables Jij. By allowing for Jij = +1 and Jij = — 1, the program accommodates 
the Ising ferromagnet and anti-ferromagnet as special cases. The multicanonical parameters 
are coded into an effective action table A{-), new spins are proposed randomly and then 
rejected with probability (1 — min[l, exp(y4' — A)]), where A' is the effective action of the 
new configuration. 

As an exercise and to check our code on exact results, we performed a multicanonical 
simulation of the 2D Ising model with < (3 < oo. We kept the time series of two million 
sweeps and measurements on a 25 x 25 lattice and verified that the finite lattice specific heat 
results of Ferdinand and Fisher are well reproduced. No difficulties are encountered 
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Figure 2: Multicanonical energy density distribution (L = 48 lattice, realization 2). 



with the multicanonical ensemble when crossing the phase transition point. To explore 
the possibility of zero temperature entropy calculations, we used Z{0) = 2^^^ as input and 
obtained = 0.61 ±0.09 for the total groundstate entropy. This corresponds to an estimate 
of 1.84±0.17 groundstates, i.e., within statistical errors we are in agreement with two. Using 
Z (0) = 2^™*^ and a time series of four miUion sweeps on a 50 x 50 lattice we obtained 2.07±0.22 
for the number of groundstates. Figure 1 is obtained from the simulation of this lattice. 

After this test we turned immediately to the 2D Edwards-Anderson spin glass. On 
lattices of size L = 4, 12, 24 and 48 we performed multicanonical simulations. Up to L = 24 
we investigated ten different realizations per lattice and, due to CPU time constraints, we 
considered only five reahzations for the L = 48 lattice. To study the simulation method, a 
rather small numbers of realizations is sufficient and to some extent desirable, as at this stage 
each single realization still deserves some individual attention. The multicanonical energy 
distribution for the second of our L = 48 realizations is depicted in Figure 2. The fall-off for 
— e < is like that of the canonical distribution at /5 = 0. For < — e < — e° an impressive 
flatness (about 800 energy entries on the lattice under consideration) is quickly achieved by 
the recursion (7). Close to the groundstate some difficulties are encountered on which we 
comment later. It should be understood that deviations from the desired constant behavior 
(6) do only influence the statistical error bars, but not the estimates themselves. Therefore, 
such deviations do not pose problems as long as they can be kept within reasonable limits 
of approximately one order of magnitude. In Figure 3 the function (3{E) is given versus E 
for the same realization as used for Figure 2. 

Tables 1-4 give an overview of our numerical results. All entries are as introduced in 
the previous sections. For /^max we take f3{E^), where it should be noted that due to our 
computational procedure P{E) is a noisy function. The final mean values and their error bars 
are obtained by combining the results from the different realizations. Different realizations 
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Figure 3: Multicanonical l3{E) function (L = 48 lattice, realization 2). 



are statistically independent and enter with equal weights. The final error bar is enlarged by 
a Student multiplicative factor, such that the probabihty content of two standard deviations 
is Gaussian (95.5% ). The dependence of the estimates on the realizations is fairly strong. 
Optimal use of CPU time is made by calculating estimates on the realizations with identical 
variance than the one obtained by the dependence on the realizations. Therefore, for the 
smaller lattices it would have been more efficient to simulate a larger number of realizations 
with less CPU time spent on each. Towards larger systems one expects self- averaging of 
energy and entropy with respect to the quenched random variables Jij, and a decrease of 
the variance between different realizations like 1/N. Our data are consistent with such a 
behavior. There is no self-averaging expected and found for quantities which relate to the 
spin glass order parameter. 

For the L = 4 lattice we also calculated the exact results by enumeration of all 2^^ 
configurations per realization. As expected, the groundstate energies reported in Table 1 
are identical with the exact values. For all other quantities of Table 1 we find reasonable 
agreement within the statistical errors. The about 2a discrepancy of our (L = 4) with 
[|] should therefore be interpreted as a statistical fiuctuation. 

In Figure 4 we plot the ergodicity time versus lattice size L on a log-log scale. The data 
are consistent with a straight line fit {Q denotes the goodness of fit [20|), which gives the 
finite size behavior 

rl ~ sleeps. (12) 

In CPU time this corresponds to a slowing down ~ \/3.2{2)^ should be remarked that a fit 
of form r£ ~ exp(cL) results in a completely unacceptable goodness of fit Q < 10~^. Still, 
the behavior (12) is by an extra volume factor worse than the close to optimal performance 
we hoped for. To understand the reasons for the rather high power, we first focused on the 
acceptance rate of the Metropolis updating, which is depicted in Figure 5 for our various 
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Table 1: L = 4 results. Each realization relies on 1,600,000 sweeps. 
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Table 2: L = 12 results. Each realization relies on 400,000 sweeps, the data points marked 
by * have four times this statistics. 
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1 


2.00 


40696 


39 


-1.3403 


0.1135 (05) 


0.18 (04) 


0.63 (11) 


2 


2.45 
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3 


-1.3958 
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3 


1.87 
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11 
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0.0789 (07) 


0.27 (06) 


0.61 (15) 


4 


2.31 


338000 


2 


-1.4028 


0.0645 (18) 


0.67 (09) 


0.99 (01) 


5 


1.94 


97587 


16 


-1.3681 


0.0851 (06) 
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6 
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10 
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23 
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5 
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9 


2.20 
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9 
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0.0868 (16) 
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10 


2.73 
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4 


-1.4028 
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Mean 
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0.081 


0.406 
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Table 3: L = 24 results. Each realization relies on 1,600,000 sweeps. 
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2.13 


425705 


2 


-1.4132 


0.0781 (5) 


0.050 (57) 


0.70 (33) 


2 


2.16 


992676 


4 


-1.4063 


0.0786 (3) 


0.043 (21) 


0.57 (38) 


3 


2.34 


2464240 


2 


-1.4080 


0.0801 (3) 


0.141 (08) 


0.94 (04) 


4 


2.33 


2004270 


2 


-1.3924 


0.0812 (8) 


0.079 (79) 


0.74 (36) 


5 


2.12 


1399679 


2 


-1.3767 


0.0914 (9) 


0.098 (06) 


0.86 (14) 


Mean 


2.22 


1457315 




-1.399 


0.082 


0.084 


0.76 


Error 


± .07 


± 516925 




± 0.010 


± .004 


± .025 


± .10 



Table 4: L = 48 results. The realizations rely on 4,000,000 — 6,400,000 sweeps. 
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Figure 4: Ergodicity times versus lattice size on a double log scale. 
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move 


4 


0.392 ± 0.026 


12 


0.196 ±0.018 


24 


0.096 ± 0.021 


48 


0.196 ±0.143 



Table 5: Metropolis "move" probability versus lattice size. 



lattice sizes. In all cases we find a drop towards about 10% in the ground state vicinity. 
This lattice size dependence is too weak to provide on its own a sufficient explanation for 
the extra volume factor encountered in (12). Therefore, we investigated next which new 
states are really accepted. We find that in the vicinity of the groundstate most accepted 
proposals do not change the energy, whereas in the disordered region the energy is changed 
in a majority of the accepted cases. Let us denote by Pmove the probability that an accepted 
new state has a different energy then the previous state. Table 5 summarizes the ground- 
state -P^o^e probabilities versus lattice size. The observed decrease is conjectured to explain 
the encountered extra volume factor. The energy random walk picture is still correct, but 
the Metropolis algorithm moves more slowly with increasing lattice size. It may be that 
rather straightforward heat bath consideration with respect to the lattice configuration as 
a whole could fix the problem, and we intend to approach this problem in future work. As 
small Pmove probabilities are limited to the groundstate neighborhood, they provide also a 
natural explanation for the difficulties which we encountered there when determining the 
multicanonical parameters. 

Figure 6 illustrates the possibility to calculate canonical expectation values for all (3 > 0. 
Energy density e = e(/3) and entropy per spin s = s(/3) for L = 48 are depicted (0 < /3 < 3). 
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Figure 5: Multicanonical Metropolis acceptance rates. 



The indicated error bars are with respect to the five different realizations. It should be 
noted that error bars at subsequent [3 values are highly correlated, because the underlying 
configurations are identical, the only difference being in the re-weighting. For small (3 the 
error bars are not visible on the scale of Figure 6 (the horizontal bar is not an error). 

We estimate the infinite volume groundstate energy and entropy from FSS fits of the 
form /l = + c/V . These fits are depicted in Figures 7. Our groundstate energy estimate 
e° = -1.394 ± 0.007 is consistent with the previous MC estimate 6° = -1.407 ± 0.008 as 
well as with the transfer matrix result [pTl e° = —1.4024 ± 0.0012. Our groundstate entropy 



estimate s° = 0.081 ±0.004 is also consistent with the MC estimate || s° = 0.071 ±0.007, but 
barely consistent with the more accurate transfer matrix result |^ = 0.0701 ±0.005. Still, 



a larger statistical sample and a careful study of systematic error sources would be needed 
to claim that there is a significant discrepancy. Our results could be improved by exploiting 
the high temperature expansion of the entropy, as it was done in . However, this would be 
against the spirit of this paper which is to explore the possibilities and limits of multicanonical 
spin glass simulations. It should be noted that the reported groundstate entropy value 
translates, even for moderately sized systems, into large numbers of distinct groundstates. 
For instance s° = 0.075 implies the following approximate groundstate numbers: 4.9 x 
10^ (L = 12), 5.8 X 10^8 (L = 24), and 1.1 x 10^^ (L = 48). 

Our results relying on groundstate data for the spin glass order parameter are of less 
satisfactory quality. Due to lack of self-averaging one has to sample many groundstates. For 
L > 12 we scaled our CPU time approximately with ~ V"^, and not ~ as would (presently) 
be required to keep the number of independent groundstates visited approximately constant. 
The number of tunneling events rir decreases with lattice size, and for L = 48 our largest 
number is only = 4 (achieved by our favorite realization). Although = 4 is a too small 
number for serious statistical conclusions, this case serves well to illuminate the problems 
as well as the achievement of the multicanonical simulation. In analogy to Figure 1, we 
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Figure 6: Energy density and entropy per spin versus /? (from the L = 48 lattices). 
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Figure 7: FSS estimate of the infinite volume groundstate energy density and entropy per 
spin. 
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Figure 8: Spin glass order parameter distribution (L = 48 lattice, realization 2). 



give in Figure 8 the spin glass order parameter distribution of this realization. Instead of 
two branches we find now four. However, there is an important difference: The number of 
tunneling events may now be small in comparison to the number of expected branches. In 
the Ising model we know that there are only two branches of the magnetization, and for 
the 50 X 50 lattice on which Figure 1 relies the number of tunneling events was Ur = 118. 
Therefore, in the Ising model we could assemble reasonable statistics for each branch. Now, 
we may have to cope with collecting a representative sample from a huge number of all 
branches. Depending on the number of tunneling events Figure 8 may look different, the 
number of groundstate "fingers" roughly proportional to the number of tunneling events. 
Figure 9 shows the groundstate state time series for the spin glass order parameter. The 
time to counts only those sweeps which emerge in a groundstate and go is the corresponding 
overlap of our independent replica. Of our four million data little more than 1% (precisely 
4,419) belong to the groundstate energy. These 1% decompose into four segments. Within 
each segment the data are highly correlated, but the segments themselves are statistically 
independent, because in-between the system tunneled each time all the way back into the 
completely disordered (3 = region. Additionally, the time series moved within each section 
occasionally out of the groundstate. This gives some amount of independence to the data 
within each segment. The numbers over each segment give the corresponding ranges in 
the total time series of our four million sweeps, and the sum of these ranges exceeds of 
course 4,419. For L = 48 the Xq ci^nd estimates suffer from insufficient statistics. Entire 
jackknife bins with go ~ exist, and error bars which may exceed the signal are implied. 
Having these limitations in mind, the smallness of our Xq data seems still to contradict the 
scaling Xq ~ [@ at the phase transition point T = 0. After improving the simulation 

method (getting rid of the constrained MC), it will be desirable to investigate a large number 
of L = 48 realizations with high statistics. 
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Figure 9: Multicanonical simulation, time series for spin glass order parameter ground state 
measurements {L = 48 lattice, realization 2). 

4 Evaluation of the Performance 

It is not entirely straightforward to compare multicanonical and standard simulations. For 
instance autocorrelation times of multicanonical simulations come out short due to the trivi- 
ality that the simulation spends most of its time at rather small effective (3 values. Therefore, 
autocorrelations are not well suited. Our ergodicity time, the average number of sweeps to 
find truly independent groundstates, is a more useful quantity. When relying on it, i.e., 
concentrating on groundstate properties alone, we should remember that this pushes the 
multicanonical simulation to its extreme. This way of calculating groundstate properties 
gives for free all properties in-between, from (3 = on. If a phase transition occurs, its study 
is also included (we noticed no particular difficulties in the 2D Ising model). For instance in 
the 3D Ising spin glass the canonical slowing down at /3c ^ is already worse than the one 
given by our equation (12). 

Although the slowing down (12) is severe, it seems to provide an important improvement 
when compared with the slowing down which canonical simulations encounter for tempera- 
tures below the bifurcation temperature. For L > 24 canonical simulations 0, ^ are unable 
to equilibrate the systems at the /?max values reported in our tables since the relaxation time 
is by far too long. Presumably due to this fact the literature focused on other questions and 
does not provide detailed FSS investigation of relaxation times. A rough estimate of the 
canonical ergodicity time may be |jl| 

^canonical ~ {^P " C) • (13) 

In this equation the scaling with L may be hidden in the L dependence of our /^max values, 
which is argued to be divergent like /3max ~ In(^), and this line of reasoning gives a slowing 
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down of the canonical algorithm like V'^^. With /3max = 2.12 (our L = 24 case) equation (13) 
leads to a canonical ergodicity time of order 10^ — 10^ when the missing constant is assumed 
to be of order one, whereas our is about 2 x 10^. 

To get a more direct handle on these problems, we performed canonical simulations 
at /3 = 1.4 for our L = 12 realizations and at /3 = 1.9 for our L = 24 realizations {/3 
values somewhat lower than the corresponding /?max averages of Tables 1 — 4 seem to be 
kind of optimal for canonical groundstate investigations). In each case two independent 
replicas were simulated and the spin glass order parameter was calculated ala (10). As our 
ergodicity time Tg and the corresponding number of tunneling events rir are nonsense for 
canonical simulations, we replace rir by a definition which allows to compare canonical and 
multicanonical simulations. Monitoring the groundstate time series q'^{t) suggests that its 
number of sign flips rij is a suitable replacement for n^. To disregard small fluctuations 
around zero, we introduce a cut-off range: A sign flip is counted when q^{t) tunnels form one 
side of the cut-off range to the other. We found that ±0.3472 gives an appropriate range, 
with L = 12 this is ±50 for q^V. According to p[ one may want to scale with L^'^ instead 
of V = L^, but the difference would not matter for our present purposes. 

Let us first report the canonical L = 12 results. We perform an identical number of 
sweeps as for the multicanonical simulations of Table 2. For each realization the canonical 
simulation finds precisely the groundstate energy reported in Table 2. In Table 6 we compare 
multicanonical and canonical flip rates. For the multicanonical case we notice that nj is 
smaller than rir reported in Table 2. This is due to the fact that subsequent, independent 
groundstate configurations agree with probability 1/2 on the sign of q^{t). For all realizations 
the canonical flip rate is smaller than the multicanonical. The actual ratio depends strongly 
on the realization. This effect can be fairly dramatic. In Figure 10 we consider realization 
^ 5 and compare its multicanonical with its canonical q^{t) time series. For the canonical 
simulation the time series never flips! Depicted are those values of q for which both replica 
are in a groundstate, and this is the case for more than a quarter of the entire time series 
(125,000 of 400,000 sweeps). In contrast to this the multicanonical simulations spent only 
8,000 sweeps, i.e., 2% of its statistics, in the groundstate, but it tunnels often (179 flips). 
Going now to replica # 6, depicted in Figure 11, the multicanonical performance is very 
similar, but the canonical simulation has now also little difficulties and flips 81 times. To 
estimate the over-all improvement factor of the multicanonical simulations is not possible 
because of the two cases with Uf = 0. If we replace these zeros by one, we obtain 46 ± 22, 
what should be considered as some kind of lower bound. 

For L = 24 we performed canonical simulations with four times the multicanonical statis- 
tics of Table 3 , i.e., 6,400,000 sweeps per replica. Nevertheless we obtain zero tunneling 
events for most of the realizations. In view of the ergodicity time estimated by our previous 
consideration, this is no surprise. By these zeros any attempt to estimate the L = 24 multi- 
canonical improvement directly is rendered hopeless: One would have to perform canonical 
simulations until at least one flip is obtained for each replica. To do so would be an obvi- 
ous waste of computer resources. The superiority of the multicanonical approach is already 
clearly established by the L = 12 results. 

When one is only interested in groundstate properties, minimization algorithms have to 



be compared. As a method simulated annealing |T8[ stands out because of its generality, 
although there are more efficient algorithms for special cases, which should be used when 
appropriate. One problem of simulated annealing is the dependence of the results on the 
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Figure 10: (a)Multicanonical versus (b)canonical simulation: Time series for spin glass order 
parameter ground state measurements (L = 12 lattice, realization 5). 
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Figure 11: (a)Multicanonical versus (b)canonical simulation: Time series for spin glass order 
parameter ground state measurements (L = 12 lattice, realization 6). 
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# 


Muca (L = 12) 


Cano (L = 12) 


Muca (L = 24) 


Cano (L = 24) 


1 


358* 


11* 


76 


7 


2 


80 


18 


4 





3 


97 


29 


20 


1 


4 


61 


2 


4 





5 


179 





44 





6 


151 


81 


4 


1 


7 


119 





52 





8 


147 


2 


4 





9 


217* 


17* 


16 





10 


483* 


112* 


8 






Table 6: Number of flips Uf (defined in the text) for multicanonical (Muca) and canonical 
(Cano) simulations on L = 12 and L = 24 lattices. The L = 12 data points marked with * 
have four times more statistics than the other L = 12 data points. 



cooling rate r = — AT/sweeps. True groundstates are only encountered for r sufficiently 



small. The r dependence of groundstate energies was investigated in p2[, and for our model 
the behavior 

e(r) = e° + cr^ with c ^ 0.5 (14) 

is indicated. To find a true groundstate, one has to reduce [e(r) — e°] to the order 1/V. 
Assuming that the constant in equation (14) is volume independent (only the lattice size 
100 X 100 was considered in [^), this translates to 

sweeps ~ \/4 ^ L^ (15) 

far worse than our equation (12). This result is kind of amazing as the multicanonical 
ensemble has eliminated directed cooling and is nevertheless more efficient. If one does not 
insist on true groundstates one can relax the condition [e(r) — e°] ~ 1/V. For instance, any 
behavior [e(r) — e°] ^0 with L oo would still give the correct density, and simulated 
annealing would slow down far less dramatic then according to (15). On the other hand, 
this would also imply a less stringent multicanonical simulation, and the final outcome of 
such a comparison is unclear. To compare the performance on lattices of fixed sizes is more 
subtle than just fixing the constant in equations (12) and (15). The reason is that the 
multicanonical simulations normally finds more than one groundstate within one tunneling 
period. This will amount to an extra factor in favor of the multicanonical simulation, which 
can be determined by calculating the integrated autocorrelation time for the groundstate 
series. Presently we did not attempt this, because at least for our L = 48 realizations the 
statistics are insufficient for such an enterprise. 

A comparison with the replica MC algorithm |^ is even less clear cut. The obtained 
estimates of the groundstate energy and entropy are in accuracy similar to ours. As one has 
to simulate many replica at many /3-values a direct comparison is impossible. Clearly, the 
results reported on slowing down are much more promising than ours for the ergodicity time. 
In particular, with our present implementation we would be unable to equilibrate a large 
128 X 128 system. On the other hand, to our knowledge the replica MC approach has never 
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been applied to the 3D Edwards- Anderson model, and one may well encounter difficulties. 
In contrast, 3D multicanonical simulations are straightforward. In fact the dimension is just 
a parameter in all our computer programs and we have already carried out various 3D test 
runs. 

Let us address the question of algorithms from a more general perspective. With a 
Metropolis type implementation our optimum performance will be bounded from below by 
a slowing down ~ V'^ in CPU time, and we are outperformed by any algorithm which can 
do better than this. For a number of important, but often highly specialized, applications 
such better algorithms exist and should be used. Cluster algorithms are presumably the 
most noticeable case of specialized high performance algorithms. The main advantage of the 
multicanonical ensemble is its generality. With this respect our method resembles simulated 
annealing , while clearly avoiding some of its disadvantages. Finally, it should be stressed 
that the multicanonical ensemble is an ensemble and not an algorithm. One may find better 
algorithms than conventional Metropolis updating to simulate this new ensemble. To try a 
combination with cluster algorithms is certainly an attractive idea. 



5 Conclusions 

Simulations of the multicanonical ensemble open new perspectives for a wide range of appli- 
cations. The present paper shows that multicanonical spin glass simulations are feasible and 
some results are very encouraging. Still, a number of technical details have remained tedious. 
Our believe is that by gaining more experience with the multicanonical ensemble more ef- 
ficient implementations of many details (like determining the parameters) will emerge. We 
hope that the present paper will stimulate further investigations. By no means we do claim 
that simulations of the multicanonical ensemble are already well understood. 

In our present implementation we find a power law slowing down ~ in CPU time. 
Quantitative comparisons with standard simulations are kind of difficult, as in awareness 
of the ergodicity problems there aims have been different from ours. In addition, a fair 
comparison has to take into account that besides groundstates our multicanonical simulation 
covers the entire range from (3 = on. This amounts to an additional factor of at least 
\/V in favor of the multicanonical ensemble. Even without doing so, the improvement 
of the slowing down of canonical simulations is remarkable. We were able to equilibrate 
the system at f3 values which are simply inaccessible to canonical simulations because of 
relaxation problems. When we constrain our attention just to the problem of identifying 
true groundstate configurations in the limit of large systems, our performance seems to be 



superior to that of simulated annealing WK Both methods share that they can address 



a very general range of applications, but a major advantage of the multicanonical ensemble 
is that the relationship to the equilibrium canonical ensemble remains exactly controlled. 
As we talk about a new ensemble (in contrast to just an algorithm), further improvements 
seems to be possible through combination with algorithms which are more efficient than 
simple Metropolis updating. In particular the possibility of exploiting cluster ideas 



could be explored. 

Qualitatively our results make clear that the multicanonical approach is certainly a rele- 
vant enrichment of the options one has with respect to spin glass simulations. The similarities 
of spin glasses to other problems with conflicting constraints [113 suggest that multicanonical 
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simulations may be of value for a wide range of investigations: optimization problems like 
the travelling salesman, neural networks, protein folding, and others. Multicanonical simu- 
lations of the 3D Edwards-Anderson model may eventually shed new light on the questions 
whether the model exhibit mean field-like behavior or some kind of droplet picture applies 
I, I, HI. In the later scenario one would expect that a single valley dominates the others 
in the thermodynamic limit, what simply means it leads down to lower energies. On the 
other hand in the mean field scenario one would expect groundstate degeneracy over many 
valleys. It seems that previous numerical work on this question is somewhat inconclusive 
as sufficiently low temperatures could not be reached without destroying thermodynamic 
equilibrium. 
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